Mangel and Bonsall Theoretical Biology and Medical Modelling 2013, 10:5 
http://www.tbiomed.eom/content/1 0/1 IS 



THEORETICAL BIOLOGY AND 
MEDICAL MODELLING 



RESEARCH 



Open Access 



Stem cell biology is population biology: 
differentiation of hematopoietic multipotent 
progenitors to common lymphoid and 
myeloid progenitors 



Marc Mangel 1 - 2 * and Michael B Bonsall- 



Abstract 

The hematopoietic stem cell (HSC) system is a demand control system, with the demand coming from the organism, 
since the products of the common myeloid and lymphoid progenitor (CMP, CLP respectively) cells are essential for 
activity and defense against disease. We show how ideas from population biology (combining population dynamics 
and evolutionary considerations) can illuminate the feedback control of the HSC system by the fully differentiated 
products, which has recently been verified experimentally. We develop models for the penultimate differentiation of 
HSC Multipotent Progenitors (MPPs) into CLP and CMP and introduce two concepts from population biology into 
stem cell biology. The first concept is the Multipotent Progenitor Commitment Response (MPCR) which is the 
probability that a multipotent progenitor cell follows a CLP route rather than a CMP route. The second concept is the 
link between the MPCR and a measure of Darwinian fitness associated with organismal performance and the levels of 
differentiated lymphoid and myeloid cells. We show that many MPCRs are consistent with homeostasis, but that they 
will lead to different dynamics of cells and signals following a wound or injury and thus have different consequences 
for Darwinian fitness. We show how coupling considerations of life history to dynamics of the HSC system and its 
products allows one to compute the selective pressures on cellular processes. We discuss ways that this framework 
can be used and extended. 

Keywords: Hematopoieitic stem cell, Multipotent progenitor, Common lymphoid progenitor, Common myeloid 
progenitor, Darwinian fitness, Natural selection, Population dynamics 



Introduction 

Hematopoiesis (the formation of blood components) is a 
highly orchestrated and dynamical process. Hematopoi- 
etic Stem Cells (HSCs) give rise, through a large array 
of successively differentiated progeny, to mature blood 
cells. While progress has been made in understanding the 
HSC system, particularly at the molecular level [1,2], Tan 
et al. (pg 82-83) [3] concluded a recent review on HSCs by 
identifying critical unanswered questions: 1) Is it possible 
to manipulate adult stem cells to increase their ability to 
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proliferate in vitro while maintaining stem cell qualities, 
so that adult stem cells can be used as a sufficient source 
of tissue for transplants and other therapeutic strategies? 
2) What are the intrinsic and extrinsic controls that keep 
stem cells from differentiating or that direct them along 
a particular differentiation pathway to form one special- 
ized cell type rather than another? 3) What are the factors 
responsible for stem cell responses to injury or damage 
that enable rapid activity and appropriate contribution 
to tissue repair and regeneration? 4) Can the 'stress' sig- 
nals' that command facultative stem cells to respond to 
tissue damage and gain specific regenerative quality be 
harnessed for therapeutic value? 

We will show that these questions can only be fully 
answered if one considers the connection between the 
needs of the organism and the HSC system, a demand 
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control system [4] in which the demand comes from the 
organism via signals to the HSC system. As noted by 
Metcalf [4] in his classic lectures, appropriate long-term 
behavior of the HSC system is essential for the health of 
the organism. Furthermore, HSCs are required to pro- 
duce the required differentiated cells without depleting 
the stem cell pool or creating damaged stem cells that lead 
to cancer. Even in the absence of oxidative damage due to 
ischemia/reperfusion, blood loss, or infections damage to 
blood cells and tissues caused by reactive oxygen species 
is a common process leading to cell apoptosis, and mature 
blood cells have a steady rate of turnover and are thus 
constantly replaced [5]. 

Thus, organismal performance is intimately connected 
to HSCs and their products [6-9]. The myeloid products of 
the HSC system are particularly important for antipreda- 
tor and foraging behaviors and general immune response 
and the lymphoid products of the HSC system are essen- 
tial for specific immune response and neonatal survival. 
Recent whole organism manipulations have demonstrated 
the veracity of this proposition (see [10] for myeloid cells; 
[11] for lymphoid cells). At the same time, when the 
organism is in a steady state (homeostasis) the HSC sys- 
tem and its products are relatively stable [12]. In addition, 
the demand is not unlimited; for example, it has been 
known for a very long time that organismal performance 
is a peaked function of hematocrit (e.g., [6,9,13]). 

In this paper, we show that the questions raised by 
Tan et al. [3] are answered most effectively if one takes 
an approach to stem cell biology based on popula- 
tion biology, which involves combining the dynamics of 
cells with considerations of fitness and natural selection 
[14]. We develop a theoretical framework that comple- 
ments existing models (Additional file 1: Table SI) to 
explore the dynamics of the HSC system. Elsewhere ([15]; 
Figure 1 here) we have developed stochastic and deter- 
ministic models of HSC cells and their associated prod- 
ucts and applied evolutionary invasion analysis [16] and 
state dependent life history theory [17-20] to show that 
understanding the dynamics of HSCs and their products 
requires asking more than whether a stem cell renews, 
symmetrically differentiates, or asymmetrically differen- 
tiates. Understanding the roles of positive and negative 
feedback is essential for predicting stem cell dynamics. By 
linking these feedback processes to stochastic population 
models (which allow uncertainties inherent in the system 
to be accounted for) we showed how well the overall mean 
dynamics of the system can be approximated a system 
of ordinary differential equations. We build on this work 
here and introduce additional concepts associated with 
population biology to the biology of stem cells, with the 
focus on the HSC system. We introduce the Multipotent 
Progenitor Commitment Response (MPCR) that charac- 
terizes the penultimate differentiation of a multipotent 



progenitor (MPP) to a Common Lymphoid Progeni- 
tor (CLP) or a Common Myeloid Progenitor (CMP), 
i.e. whether they follow a myeloid or lymphoid track. 
Although we recognize that within the myeloid track there 
is another decision towards a granulocyte-macrophage 
progenitor or megakaryocyte-erythrocyte progenitor (see 
Additional file 1). We show how the fitness (survival and 
reproduction) of the organism shapes the MPCR, thus 
providing an approach for modeling the demand control 
nature of the HSC system. 

Methods 

We begin first by describing, in summary here, with 
details in Additional file 1, the dynamics of the stem cells 
and their descendants, after which we describe the com- 
ponents of fitness (survival and the reproduction) and 
their dynamics. We then couple the two together. 

Dynamics of stem cells and their descendants 

In the Additional file 1, we derive the computational 
model given below and in Figure 1 we provide a graphical 
representation of the full mathematical model. We thus 
consider stem cells, with concentration denoted by [5], 
Multipotent Progenitor cells, denoted by [MPP], Com- 
mon Lymphoid Progenitors, denoted by [ CLP], Common 
Myeloid Progenitors, denoted by [ CMP], fully differenti- 
ated lymphoid cells, denoted by [L] (measured in num- 
bers per milliliter), and fully differentiated myeloid cells, 
denoted by [M] (measured in numbers per nanoliter). 
Central to these dynamics we assume that the stem cell 
niche can support at mostJC stem cells and that in absence 
of all other feedback (described below), the dynamics in 
the niche follow Gompertzian kinetics (justified in [15]). 

The dynamics are described by the following set of 
coupled ordinary differential equation: 

d[S] 

[S] -log(K/[S] )(r s - r p ,® p ,([L] , [M] )) 



dt 



d[MPP] 
dt 



d[ CLP] 
dt 



d[ CMP] 
dt 



d[L] 

dt 
d[M] 

dt 



x <J> s ([L],[M])- Ms [S] 



(1) 



= [S] -log(K/[S] )(r, + 2r p ^([L] , [M] )) (2) 

x 3> S ([L] , [M] ) + (k - r dMPP )%([L] , [M] ) 
x[MPP] -fi p [MPP] 

= r dMPP <S> p (L,M)p([L] , [M] )Q N [MPP] 
-r C Lp[CLP]-ii C Lp[CLP] (3) 

= r dMPP <S> p {L,M)(\ - p([L] , [M)] )Q N [MPP] 
- r CM p I CMP] -ficMpi CMP] (4) 

= r C Lp[ CLP] +{ri - in - ^hIv(t)>v th )[L] (5) 

= r cmp [ CMP] +{r m - ii m )[M] . (6) 
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Figure 1 A diagrammatic derivation of Eqns 1 to 6 (details given in Additional file 1 ). a) In the most general case, we consider stem cells (S), a 
series of Multipotent Progenitor Cells (MPP), a Common Lymphoid Progenitor (CLP) and a Common Myeloid Progenitor (CMP). CLPs give rise to B, 
NK, and T cells; CMPs give rise to Erythrocytes (E), Granulocytes (G), and Platelets (P). We denote the total numbers of lymphoid and myeloid cells by 
Land M respectively, rates of differentiation by r. (with subscript indicating the cell type involved), rates of development of MPP cells by k., feedback 
from fully differentiated cells on those rates by and rates of cell death by /x..The feedback functions have the property that they are 1 when 
stem cell or fully differentiated cell numbers are low and decline as stem cells or fully differentiated cells increase. Thus, for example, stem cells 
renew (one stem cell becomes two) at rate r s 4> s (/, m) when the concentrations of lymphoid and myeloid cells are / and m respectively, 
asymmetrically differentiate (one stem cell becomes two stage-0 progenitors) at rate 2rrf$pi(J,m)& s (l,m), symmetrically differentiate (one stem cell 
becomes a stem cell and a stage-0 progenitor) at raterp<J> s (/, m), and die at rate /x s - Similar interpretations hold for other transitions. The Multipotent 
Commitment Response (MPCR), denoted by p(/,m), is the probability that a MPP in its final stage commits to the lymphoid route, b) To focus on the 
MPCR, we combine all of the fully differentiated cells into lymphoid and myeloid classes (L and M) and use Michaelis-Menten-like arguments to 
compress the MPP class into a single stage, assuming that steady states of intermediate stages are rapidly reached, characterized by combination of 
rate constants £2/v. 



where <&,([£], [AT] ), 4y ([L] , [M] ) and %([L],[M]) 
are feedback functions (see below) for the activity of stem 
cells (s), the asymmetric differentiation of stem cells (p'), 
and the activity of MPP cells (p), respectively. p([ L] , [M] ) 
is the demand control function (see below) that describes 
the probability of an MPP cell differentiating into a lym- 
phoid or myeloid progenitor. Q.}q[MPF\ is the survival 
of MPP cells from initial differentiation through to ulti- 
mate differentiation into CLP or CMP cells (see Additional 
file 1 for its full derivation), /i^ characterizes the addi- 
tional mortality when the immune system is activated 
and I a> b is an indicator function that is 1 if a > b 
and 0 otherwise. All other parameters are defined in 
Table 1. 

Feedback control, which requires nonlinear dynamics, 
is essential for the growth and regeneration of tissues. 
A recent model of genetic products [21] to character- 
ize the erythroid-myeloid lineage decision, shows how 
nonlinearities arise, de Graaf et al. [22] showed that the 



platelet concentration can regulate HSCs, through the 
concentration of Thrombopoietin (TPO). Other authors 
have shown that Lkbl, a kinase enzyme best known as 
a tumor suppressor, provides feedback control of HSCs 
[23-26]. Such work provides the empirical basis for our 
modeling. 

To incorporate feedbacks we extend Lander et al.'s [27] 
approach: if r denotes a generic reaction rate constant 
and [/] the concentration of fully differentiated cells, they 
model the reaction rate with feedback control is x+TTJTJ- 
We have adapted this framework for the HSC system, with 
changes. First, two kinds of differentiated cells provide 
feedback; we use / and m to denote the total concentration 
of lymphoid and myeloid cells. Second, there is poten- 
tially different feedback on the activity of stem cells, the 
asymmetric differentiation of stem cells, and the activity of 
MPP cells. Third, the system needs to be active when there 
is a shortage of either lymphoid or myeloid cells. Thus 
we set: 
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Table 1 Variables, parameters, their interpretation, and values 



Symbol Interpretation Value 

t Non-dimensional time 1-3500 

[5] Concentration of stem cells at timer Eqn 1 

[MPP] Concentration of Multipotent Progenitor (MPP) cells at time f Eqn 2 

[CLP] Concentration of Common Lymphoid Progenitor (CLP) cells at time f Eqn 3 

[CMP] Concentration of Common Myeloid Progenitor (CMP) cells at time f Eqn 4 

[L] Concentration of fully differentiated Lymphoid (L) cells at time f Eqn 5 

[ M] Concentration of fully differentiated Myeloid (M) cells at time r Eqn 6 

K Maximum number of stem cells in a niche 10 

r s Maximum rate of stem cell self-renewal 2.5 

fp< Maximum rate of stem cell asymmetrical division 0.001 

<&'p(lL] , [M]) Feedback control from fully differentiated cells to asymmetric division Eqn 9 

*s(U] , [Afl ) Feedback control from fully differentiated cells to stem cell self-reneval Eqn 7 

* P (U] , [M] ) Feedback control from fully differentiated cells to symmetric division Eqn 8 

/x s Rate of stem cell death 0.004 

X Rate of MPP multiplication 0.25 

Hp Rate of MPP cell death 0.02 

£2/v Combination of intermediate multipotent progenitor rate constants 1 .0 

rap Rate of division of CLP into fully differentiated lymphoid cells 0.01 

/lap Rate of CLP cell death 0.001 

n Rate of multiplication of lymphoid cells 0.025 

IM Rate of lymphoid cell death when immune system is not activated 0.028 

/x/ t Additional rate of lymphoid cell death when immune system is activated 0.01 

/ 0> b Indicator function for the inequality =1 if a > b, 0 otherwise 

v th Threshold concentration for pathogens to activate the immune system 0.025 

r C MP Rate of division of CMP into fully differentiated myeloid cells 0.01 

tiCMP Rate of CMP cell death 0.001 

r m Rateof multiplication of myeloid cells 0.0 

/x m Rateof myeloid cell death 0.01 

/ Value of [L] varies 

m Value of [M] varies 

</> s /(/) Feedback control of fully differentiated lymphoid cells on stem cell activity Eqn 1 0 

0 sm (m) Feedback control of fully differentiated myeloid cells on stem cell activity SimilartoEqn 10 

<fr p l(l) Feedback control of fully differentiated lymphoid cells on symmetric renewal Eqn 1 1 

&pm(m) Feedback control of fully differentiated myeloid cells on symmetric renewal Similar to Eqn 10 

ipp'l(l) Feedback control of fully differentiated lymphoid cells on asymmetric renewal Eqn 1 1 

0p'm(m) Feedback control of fully differentiated myeloid cells on asymmetric renewal Similar to Eqn 10 

a,j Feedback parameter in <p S j(l) 1 0 

a p i Feedback parameter in 0 p /(/) 100 

a p 'i Feedback parameter in tp p 'i(l) 20 

a !m Feedback parameter in 0 sm (m) 0.1 

a pm Feedback parameter in 0 pm (m) 0.001 

a p i rn Feedback parameter in tp p r m (m) 0.2 

a Coefficient in MPP Commitment Response (MPCR) Varies 
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Table 1 Variables, parameters, their interpretation, and values (Continued) 


Symbol 


Interpretation 


Value 


Y 


Exponent in MPCR 


Varies 


K\ 


Density of lymphoid cells in homeostasis 


30 


Km 


Density of myeloid cells in homeostasis 


30,000 


Ph 


Fraction of lymphoid cells in homeostasis 


Eqn 15 


Af([M] 


Rate of accumulation of fitness when myeloid cell concentration is [M] 


Eqn 17 




Fitness accumulated to time f 


Eqn 18 


e 


Ratio of organismal to cellular time scale 


0.05 


5(0 


Survival to time r 


Eqn 21 




Total rate of mortality when myeloid cell concentration is [M] 


Eqn 20 


MeO 


Myeloid independent rate of mortality 


0.05 


Mel 


Myeloid dependent rate of mortality 


5.0 


lx,{[V\) 


Additional rate of mortality when concentration of infectious agents is [ V] 


Eqn 20 


ina 


Coefficient of [ V] in additional mortality 


0.02 


my 


Coefficient of [ V] 2 in additional mortality 


0.002 


r v 


Replication rate of infectious agents 


0.05 


Q 


Clearance rate of infectious agents by lymphoid cells 


0.05 


Cm 


Clearance rate of infectious agents by myeloid cells 


0 


vo 


Concentration of infectious agents at the start of an infection 


1 


Vth 


Concentration of infectious agents below which additional lymphoid mortality does not occur 


0.025 



(Parameter values are a canoncial fixed set, arbitrarily chosen, to illustrate the general principles of an MPCR). 



4> s (/, m) = max[ 4> s i(l), 4> sm (m)] 
<t> p (l, m) = max[ (j> p i{l), (j> pm (m)] 
<iy (/, m) = max[ 4>p>i(l), (p p > m (m)] 



(7) 
(8) 
(9) 



where </> s /(0) = </> sm (0) = 1 (representing the feedback of 
lymphoid and myeloid cells on stem cell activity), etc, and 
all <t>ij are decreasing functions of their arguments, as in 



</>d(D 



cppliD 



l + a sr l 
1 



1 + a p i ■ I 
1 + awi ■ I 



(10) 
(11) 
(12) 



where a s i, a p i, a p q are parameters. A similar form is used 
for the feedback control from myeloid cells. This is the 
simplest form of feedback between the whole organism 
and the bone marrow stem cell system; see [28] for alter- 
natives. 

In a demand control system the probability of a MPP 
cell ultimately differentiating into a CLP or CMP cell must 
depend upon the state of the organism. That is, the cur- 
rent densities of myeloid and lymphoid cells determine 
the appropriate response. We choose a functional form 



that is widely used in population biology and similar to 
Michaelis-Menten enzyme kinetics; 

~(f) y 



p(l,m) = 



(13) 



! + «(?)" 

where p(l, m) represents the demand control function and 
a and y are parameters that describe this specific asymp- 
totic Michaelis-Menten function. We let K m : ki denote 
the ratio of myeloid to lymphoid cells in homoeostasis. 
If Ph denotes the value of p {I, m) in homeostasis then on 
average we have: 

Ph = (14) 

and from Eqn 13 we have 

a(K m /Ki)Y 

Ph = — — (15) 

1 +a(K m /Ki)y 

We view the unknowns in this equation as the two param- 
eters a and y from which we find 



Ph 



Ph 



(K m /KlY 



(16) 



Eqn 16 determines a curve in the y-a plane and every 
value on this curve provides the same value of p (I, m) = 
Ph- However, when out of homeostasis, the value of a(y) 
that determines Eqn 16 has a profound effect on the 
fate of HSC descendants. In Figure 2a we show the rela- 
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tionship between the parameters a and y of the MPP 
commitment response when homeostasis corresponds to 
1 lymphoid cell per 1000 myeloid cells, a typical ratio for 
humans. Each point on the lines in these panels corre- 
spond to a particular value of the pair (y,a) consistent 
with the number of cells at homeostasis. In Figure 2b, 
we show how p(l,m) varies as m/l varies for three val- 
ues of y. As y increases, the MPP commitment response 
becomes more sensitive to variation in the densities of 
myeloid and lymphoid cells. These three curves represent 
just a few possibilities in the infinite space of functional 
responses corresponding to Eqn 13 and raises the ques- 
tion: how can we predict which response an organism will 
use (i.e., where on the curve in Figure 2a will a population 
of organisms sit)? We discuss this below. 

How an organism goes out of homeostasis depends 
upon its environment. For example, in an environment 
when wounds occur frequently, we anticipate the m/l will 
be lower than the value in homeostasis while in an envi- 
ronment when infection occurs frequently we anticipate 
that m/l will be greater than that value in homeostasis. 

In this paper, we are interested in p(l, m), and in partic- 
ular how natural selection affects it, in light of the envi- 
ronment of the organism. To do this, we need to couple 
Eqns 1-13 to organismal fitness, which is what separates 
our work from all that has come before it. We consider it 
advisable to step back from modeling a particular situa- 
tion, but instead consider the more general properties that 
connect the needs of the organism with the activity of the 
bone marrow stem cell system. 

The components of fitness and their dynamics 

The representation of genes in subsequent generations is 
determined by survival and successful reproduction of the 



focal organism. Regarding the latter, we assume that the 
rate at which successful reproduction occurs {Af ([M] )) 
is a function of myeloid cells, justified by the long recogni- 
tion that organismal performance is a peaked function of 
hematocrit. Following figure two in [9], we set 



Af([M] ) = a 0 +ai-M + a 2 -M 2 



(17) 



provided this expression is positive; otherwise we set 
A/([Af]) = 0, where fl 0 = -1.034565, a\ = 0.001527, 
and fl 2 = -0.0000002864. The peak of Af([M] ) occurs 
at [M]* = 2666 (Figure 3). Hematocrit is the fraction of 
myeloid cells in blood; if too great then fitness is impaired 
as an organism has few lymphoid (immune) cells. If too 
low then there are insufficient myeloid cells to support 
the oxygen needs of the organism. Physiologically, the 
parabolic hematocrit function emerges through the rela- 
tionship between the oxygen carrying capacity of blood 
and blood flow. When myeloid cells are low (M — > 0), 
blood flow is maximal but oxygen concentration is low. In 
contrast as M — > oo, oxygen concentration is maximal but 
blood is too viscous to flow. Hence, this proximate phys- 
iological explanation and consequently the ultimate evo- 
lutionary fitnees cost lead to a peaked (parabolic) hemat- 
ocrit function. The numerical relationship between repro- 
ductive rate and myeloid cell concentration described in 
equation (17) captures this parabolic shape. 

We let J-(t) denote lifetime fitness accumulated to time 
t, S(t) denote survival to time t, and e << 1 a scaling 
parameter that relates the organismal and cellular time 
scales. Then 



= e ■ S(t) ■ Af([M(t)]) 

at 



(18) 




-2.5- 



e -3 - 



7=1.5 




500 1000 1500 2000 2500 3000 3500 4000 



Y m/l 

Figure 2 a) The relationship between the parameters a and y of the stem cell commitment response when homeostasis corresponds to 1 
lymphoid cell per 1000 myeloid cells, b) Different values of y affect how the MPCR varies with changes in the number of lymphoid and myeloid 
cells. In the presence of high numbers of myeloid cells, the demand response is to drive the MPPs to make more lymphoid cells. 
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Figure 3 We assume that the rate at which successful reproduction accumulates, Af(m) is a parabolic function of the density of myeloid 


cells m. 





To determine survival, we assume that uninfected indi- 
viduals have a per unit time rate of mortality with myeloid- 
independent and myeloid-dependent components so that 
the total rate of mortality is (Figure 4) 



fi e ([M] ) = /x e o - 



Mel 

[M\ 



(19) 



Although we focus on non-fatal diseases here, such dis- 
eases can still increase mortality rate, e.g. by reducing 
the effectiveness of flight responses. Hence we assume 
that the additional mortality induced by the pathogen 
is 



m([ v\) = wo[v] +mi[vY 



.that 



i- = - e . (^([M]) + ^[V] )S 
at 



(20) 



(21) 



Finally, we incorporate the dynamics of infectious 
agents. We assume that in the absence of immune 
response the growth of the infectious agent is exponential 
with rate r v and that lymphoid and myeloid cells clear the 
infection at rate c; and c m respectively. In this model, again 



for simplicity, we ignore memory in the immune system. 
Thus, if V(t) denotes the density of pathogens 



dV 

— = [ r v - Cl L(t) - c m M(t)] -V(t) 
at 



(22) 



Integration of Eqns 17-22 forward in time, conditioned 
on p([L] , [M] ) and a sequence of wound and/or infec- 
tion events, allows us to compute the lifetime reproductive 
success associated with that particular MPP commitment 
response. 

Eqns 1-22 are a set of deterministic ordinary differ- 
ential equations that link the behavior of the stem cell 
system with the needs of the organism. However, organ- 
isms in nature experience wounding and infection in a 
quasi-random manner. We account for this in the follow- 
ing way. Imagine that there are K w and Ki times at which 
wounds or infections can occur (these values could, of 
course, be random variables but we treat them as fixed in 
this paper, only for purposes of simplicity) and then deter- 
mine a sequence of times T w (k w ), k w = 1,2, ....,K W and 
Ti(ki),ki = 1, 2, ....,Ki at which either a wound or infec- 
tion occurs (in principle both could occur at one time). To 
illustrate the ideas, we assume that when a wound occurs, 
myeloid cells drop by 40% and that when an infection 
occurs, the infectious agent increases to the level vq. These 
occur instantaneously and we then continue with the solu- 
tion of the differential equations. For the results shown 
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here, we assume that K w = Ki = 7 and that the times are 
uniformly distributed over the interval between day 0 and 
day 1500. 

Results 

As introduced above (Figure 2a) a homeostatic p/, ratio 
of fully differentiated myeloid to lymphoid cells specifies 
a curve in the y — a plane in which all points on this 
curve are consistent with pi, but as illustrated in Figure 2b, 
different values of y (and thus a) will lead to different 
MPP commitment responses when the organism is out of 
homeostasis. Thus, to understand how natural selection 
will act on p(l, m) requires dealing explicitly with compo- 
nents of fitness and linking them to the dynamics of the 
HSC system. 

To explore this hypothesis, we assume that the rate 
at which the organism accumulates fitness, Af(m) is a 
parabolic function of myeloid cells (Figure 3) and that the 
rate of mortality is a declining function of myeloid cells 
(Figure 4) and an increasing function of the density of 
infectious agents (Eqn 22). 

In a 'deterministic' or laboratory environment with nei- 
ther wounding nor infection, organisms still die, so that 
survival declines with age (Figure 5a) and fitness accu- 
mulates but ultimately saturates because of the declining 
survival (Figure 5b). The strength of selection on p(l, m) 



will depend upon how fitness varies with y. We show this 
in Figure 6 fitness as a function of y for the determin- 
istic case, wounding only, infections only, and wounding 
and infections. Clearly there is little selection on y in the 
laboratory case or of wounding only. Selection does occur 
when there are infections, with larger values of y, leading 
to a more responsive MPCR. 

The alternative to varying y and holding the environ- 
ment at one stochastic realization is to hold y constant 
and consider multiple realizations of the stochastic envi- 
ronment. We show the results of such an approach in 
Figure 7, for y = 1. Using these, we can compute both a 
mean and variance for the selection acting on y. 

Discussion 

Here we have developed and analysed a theoretical frame- 
work for linking the population biology of the hematopoi- 
etic stem cells to the demands of the individual. We have 
introduced the notion of the MPCR (Multipotent Pro- 
genitor Commitment Response) as the response that the 
describes the penultimate decision of stem cells before 
commitment to either a myeloid or lymphoid lineage. We 
use this response to investigate the control dynamics of 
a hematopoietic stem cell system and show that different 
values of the 'shape' parameters that describe the MPCR 
give a range of optimal response. Below we discuss the 
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Figure 5 Even in a laboratory environment, without wounding or infection, organism do not live forever, so that survival declines with 
age (panel a) with the consequence that accumulated fitness saturates. 



implications of this on the evolutionary dynamics of the 
HSC system and, more broadly, for developing a theory of 
stem cell systems based in population biology. 

The meaning of a flat fitness surface 

The first derivative, |£, of the accumulated fitness with 

dy 

respect to y is a fundamental measure of the strength of 



natural selection on y. As the derivative becomes smaller, 
the strength of natural selection becomes weaker, with 
the implication that a wider spectrum of values of y will 
provide equal values of fitness. But, as we have shown 
in Figure 2b, differing values of y will lead to consid- 
erably different MPCRs, and thus kinetics of the HSC 
descendants following an external challenge such as a 
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Figure 6 Plotting lifetime accumulated fitness as a function of y allows us to understand the strength of selection on y as determined by 
the environment in which the organism lives. 
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transplant or perturbation, with the prediction that if one 
uses animals with little evolutionary history of wounding 
or infection, a wide range of HSC dynamical responses 
is expected. For instance, among 44 laboratory mice, 
Abkowitz et al. [29] observed seven different patterns of 
donor cell dynamics following hematopoeitic stem cell 
transplant experiments, suggesting that there is individ- 
ual heterogeneity in the parameters of the MPCR, as we 
would predict. 

In previous work [15] we showed that the differen- 
tial equations used here are a good approximation for 
the mean of underlying stochastic system. Understand- 
ing the limitations imposed by stochastic fluctuations on 
the feedback in our model [30] is an important next 
step because the comparisons of models and data will 
require a framework to account for process stochasticity 
and observation error. Developing appropriate state-space 
approaches to understand the consequences of these het- 
erogeneities and nonlinearities on stem cell dynamics is 
clear and obvious future step. 

Recently, Huang [31] proposed that a systems biology 
of stem cells should consider not only the link between 
observation and process through state-space approaches 
(e.g. [32]) but also the consequences of nonlinearity and 
non-genetic heterogeneity in cell systems. Here, to start 
to understand how nonlinearity effects influence stem cell 
behavior we have introduced the notion of a MPCR under 



the control from the needs of the organism. Furthermore, 
under a paradigm of stem cell heterogeneity, we should 
anticipate variation at multiple levels within any stem cell 
system; particularly so with individual variation in the 
MPCR (y and a). 

How this model can be extended 

There are a number of extensions that go beyond the 
current work for linking population biology and stem 
cell systems. For instance, experimental validation of the 
MPCR would require repeated cell count measures of 
long-term HSCs, short-term HSCs and a range of HSC- 
derived products. Given such experimental data we envis- 
age that it would be possible to assess goodness of fit 
between an MPCR model and data (and also character- 
ize unexplained heterogeneties) using computational and 
statistical methods (Bonsall and Mangel, unpublished). 

Furthermore, our framework could be extended to study 
the consequences of transplant or perturbation effects A 
perturbation experiment can be modeled by starting the 
HSC system in its steady state and then reducing the num- 
ber of lymphoid or myeloid cells and then integrating 
Eqns 1-22 forward. In addition to predicting the kinetics 
of fully differentiated cells, we can predict the activity of 
the stem cells following the perturbation. As described in 
the Additional file 1, a transplant experiment involves the 
addition of the dynamics of host cells, which are dying, 
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and allows us to predict the fraction of fully developed 
lymphoid or myeloid cells. Our approach can also be used, 
in conjunction with evolutionary invasion analysis as in 
[15], to predict the outcome of a partial transplant of stem 
cells, in which stem cells with different kinetic properties 
are transplanted into an organism with a healthy but not 
vibrant stem cell system and for which we wish to pre- 
dict whether the transplanted stem cells will co-exist with 
the host stem cells, overtake them, or be extinguished by 
them. 



Computing p(l,m) from more fundamental principles: the 
fitness control hypothesis 



not infected at time s and that [L(s)] = I and M[s] = m. 
Then 



We began with the assumption that p(l, m) 



but there are an infinite number of functions that would 
have the properties we seek: bounded by 1 and increasing 
as the the ratio m/l increases. Given experimental data, it 
is possible to use Bayesian methods to determine p(l, m) 
(Bonsall and Mangel, unpublished), but we may also ask 
if p{l,m) could be constructed from some more funda- 
mental principles. One way to do approach this woul be 
through state dependent life history, as implemented by 
stochastic dynamic programming (SDP) [17-20], which 
we explain now but leave the details for a subsequent 
paper. 

Our results point to organismal fitness being a func- 
tion of fully differentiated lymphoid and myeloid cells, 
T{l,m,t), which we have modeled (Eqns 17-22). Begin- 
ning there, consider the question from the perspective 
of the organism, rather than from the perspective of the 
dynamics of cells. Suppose that s denotes the time scale 
for the organism, so that in the interval s to s + 1 the 
time At > > 1 elapses in the cellular model. During that 
elapsed time, we assume that up to V fully differentiated 
cells are produced and the challenge for the organism is 
to allocate these fully differentiated cells to lymphoid or 
myeloid cells. That is, if SI and Sm denote the number of 
fully differentiated cells produced during At, subject to 
Sl+Sm < V, then we seek the value of V and combination 
of &i and S m that maximizes organismal fitness. 

To do this, we let F(l, m, s) denote the accumulated life- 
time reproduction from time s until time S given that 
[L(s)] = I, [M(s)] = m and the organism is not infected. 
Similarly, we let Ft{l,m,v,s) denote accumulated life- 
time reproduction under the same conditions about lym- 
phoid and myeloid cells and that the density of infectious 
agents is [ V(s)] = v. Since fitness cannot be accumu- 
lated after time S, we have the end condition F(l, m, S) = 
Fi(l, m, v, S) = 0. Methods of SDP, along with value iter- 
ation, allow us to compute the stationary values of these 
functions fors << S [15]. 

Such a method allows us to represent the optimal pro- 
duction of CLP and CMP cells given that the organism is 



V*(l, m, s) = SqI* (I, m, s) + Som* (I, m, s) 



(23) 



is the optimal production of MPP cells at time s, given 
the current state of the organism. Similar arguments gen- 
erate the optimal production of MPP cells, distributed as 
CLP or CMP for the case of the organism being infected at 
time 5. 

In the stationary state, fitness is only a function of /, m 
and v; we denote these values by F{l,m) and Fi{l,m). 
Similarly, the analog of the quantities in Eqn 23 are 
8ol*(l,m), &om*(l,m) and V (l,m). Thus, for an organ- 
ism that is currently not infected, the fraction of CLP cells 
produced is 



p(l,m) = = 



S 0 l* (I, m) 



SqI* (I, m) + Som* (I, m) 



(24) 



This equation provides an explicit form for p {I, m) 
based directly on fitness of the organism, with a similar 
one holding for the situation in which the organism is 
infected. 

Connection to empirical studies 

Our work complements the rapid recent development 
of understanding how gene products that regulate HSCs 
operate [33]. However, translating that understanding into 
a functional perspective at the organismal level is less 
well developed. Recent work suggests that epigenetic fac- 
tors [34] are important in controlling HSC behavior [35], 
but the scale of those factors is unknown, and evidence 
suggests that HSC behavior is regulated by circadian oscil- 
lations modulated by photic cues [36], as is the behavior 
of whole organisms [37]. Other the hand, it is less clear 
how the differentiated T and B cells are regulated in adap- 
tive immunity [38]. DNA damage increases in HSCs with 
age and such damage is more poorly repaired in older 
individuals [39,40]; lymphocytes also age [41]. Rossi et al. 
[42] argue that the reduction in adaptive immune capacity 
with age and the increased incidence of myeloproliferative 
diseases both have origins in changes in the HSC sys- 
tem, presumably due to increased damage of HSCs that 
could include telomere shortening or other genome main- 
tenance failures, and other damage to DNA or proteins 
[43-46]. 

Thus, currently much is known about the mechanism 
of HSC system and its descendants [47,48] but much less 
is known about how natural selection has shaped these 
mechanisms, even though stem cells are units of nat- 
ural selection [49]. Using classic ecological methods of 
allometric analysis (e.g., [50]) and compartmental mod- 
eling, it has recently been argued [51-53] that the fun- 
damental architecture of the HSC system is unchanged 
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across mammalian species (also see [54]), which sug- 
gests a long evolutionary history determining the role 
of HSCs in organismal ecology. Recognizing that evo- 
lutionary considerations are essential for understanding 
every problem in biology [55], our goal here has been 
to understand how the penultimate route of differentia- 
tion towards a common myeloid or lymphoid progenitor 
cell [56-58] can be characterized both dynamically and 
evolutionarily. 

One common way for approaching an understanding 
of the dynamics of HSCs and their products is through 
irradiation and transplant experiments (e.g., [59]), which 
also have important implications for bone marrow trans- 
plants [60,61]. The typical experiment proceeds as follows 
[61]: a host animal receives a dose of irradiation sufficient 
to destroy the blood system (about 900 rad) but which 
leaves the rest of the organ systems intact (the intestinal 
system requires 1200-12,000 rad for irreversible damage, 
and the brain more than 12,000). This dose will kill the 
animal within 2-3 weeks, which gives a sense of the rate 
of decline of HSC descendants. Cells of a donor are iso- 
lated from the bone marrow by flushing sterile medium 
through the bone marrow and then injected into the host 
animal. The donor cells first circulate in the blood stream 
and then repopulate the bone marrow. The success of 
such transplants, from the perspective of the whole organ- 
ism is demonstrated when a skin graft from the donor 
to the host is not rejected. In the period following trans- 
plantation, the donor HSCs produce progenitors that then 
differentiate into myeloid or lymphoid cells. Each differen- 
tiation that leads to a commitment of one route or another 
involves complex signaling [56,57,62,63] and the signals 
themselves must be shaped by natural selection. We call 
this a transplant experiment. 

A less commonly used approach, but equally impor- 
tant, is a perturbation experiment in which an animal 
is challenged in a way that reduces its complement of 
erythrocytes, platelets or granulocytes (commom myeloid 
progenitor (CMP) descendants) or lymphocytes (com- 
mon lymphoid progenitor (CLP) descendants) and HSC 
activity is observed subsequent to the perturbation. For 
example, Cheshier et al. [10] bled mice on days 3, 6, and 
9 and then sacrificed them on day 10 and measured the 
markers of HSC activity (see [64] for a similar experiment 
regarding the epidermis). Baldridge et al. [11] showed that 
quiescent HSCs were activated in response to an infec- 
tion, another perturbation experiment. The models that 
we have described here can be used to predict dynam- 
ics of stem cells and their descendants for both kinds of 
experiments. 

Conclusions 

The use of quantitative models to understand the HSC 
system can be traced to the classic work of Till et al. [65], 



who used branching processes to interpret their experi- 
mental results on the variation of spleen colonies formed 
after a transplant experiment but evolutionary consider- 
ations rarely appear in models in stem cell biology [15]. 
The concepts developed for understanding how popula- 
tions of individuals respond to these abiotic and biotic 
processes have colloraries at the cellular level. In partic- 
ular, the size of stem cell systems is determined by the 
availability of molecular resources, feedbacks from differ- 
entiated cells, the effects of the epigenetic environment 
and the size of the environment capable of supporting 
stem cells [66]. Finally, population biology forces us to 
recognize the distinction between typological thinking 
(that all individuals of the same species are identical with 
constant characteristics and responses) and population 
thinking (that variation is real because every individual is 
unique and individual variation is central in ontogenetic 
and evolutionary history) [67]. While there are verbal and 
pictorial models (e.g., [5,68,69]) of how signaling between 
differentiated cells and HSC may work, we have provided 
the first general quantitative theory linking how feedbacks 
through organismal need shapes the performance of the 
HSC system (cf [70]). 

In conclusion, our work raises the critical question 
of how we connect the MPCR with the vast under- 
standing on how signalling shapes HSC products. For 
instance, appreciating how the MPCR links to T-cell 
specification [71] and drives differential levels of the sev- 
eral well characterised precursors (of which lymphoid 
primed MPPs are one) requires further consideration of 
the both the key qualitative and quantitative positive and 
negative feedbacks in these sorts of signalling pathways. 
Similarly, linking the ideas of demand feedback to ery- 
thropoietin production [72], that differs between foetus 
and adult, requires systematic approach to combining 
the dynamical drivers with the lifetime fitness of the 
organism. 

While the broad evolutionary ecological of HSC activ- 
ity such as mounting an immune response to infection 
are well-known to have a cost on reproductive fitness 
(e.g., [73]), appreciating the ecological setting suggests 
that standard laboratory animals (bred for many gener- 
ations under refined conditions) may not be the most 
appropriate organisms for the sorts of empirical tests 
of the theory we propose (also see [74] on metabolic 
morbidity of laboratory rodents). Studies on the evo- 
lutionary ecology of model organisms (see [75] for an 
approach using zebrafish) should have greater prominence 
in the understanding the population biology of stem cell 
responses. 

Although raised almost two decades ago [76] this sys- 
tems biology approach to stem cell dynamics must include 
evolutionary thinking. Approaching stem cell biology as 
population biology has much to offer both fields. 
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